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1. 


INTRODUCTION 


Sparse  matrix  methods  have  been  important  in  practical  linear  programming 
from  the  very  earliest  implementations  of  the  revised  simplex  method. 

These  methods  have  become  even  more  important  as  the  size  of  problems 
presented  for  solution  has  continued  to  grow  until,  at  present,  linear 
programs  with  a  few  thousand  constraints  are  quite  commonplace.  Fortun¬ 
ately,  and  not  surprisingly,  the  densities  of  the  constraint  matrices 
decline  with  the  increase  in  dimension  thus  keeping  the  amount  of  data 
storage  and  arithmetic  within  possible  proportions. 

Techniques  for  handling  sparse  matrices  are  important  in  two  major  areas 
in  linear  programming:- 

1.  Storage  and  retrieval  of  the  problem  data.  There  are  several 
techniques  for  storing  sparse  matrices  in  use  at  present  (see  Smith  Ci^l 
and  de  Buchot  C63). 

2.  Maintaining  and  applying  the  basis  inverse  or  substitute  Inverse. 

This  paper  is  concerned  with  the  second  area.  Readers  will  be  assumed  to 
be  familiar  with  the  simplex  method  using  the  product  form  of  inverse 
(see  Dantzig  C7j)  and  the  Gaussian  elimination  for  solving  systems  of 
equations  (see  Forsythe  and  Moler  £103). 

In  recent  years  it  has  been  realized  that  the  sparsity  of  the  inverse  repre¬ 
sentation  produced  by  the  linear  progransning  inversion  routine  can  be 
drastically  improved  by  using  the  Gaussian  or  Elimination  form  of  inverse 
(E.F.I.)  rather  than  using  the  Gauss-Jordan  or  product  form  of  inverse 
(P.F. 1.).  This  scheme  was  first  advocated  by  Markowitz  Cl  13  and  subsequently 
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for  tho  spocial  caso  of  staircase  matrices  by  Dantzig  VQ.  Later  still 
Dantzi^;,  Harvey,  McKniRht  and  Smith  C93  oxporimon tally  demonstrated  the 
superiority  of  the  L.r.I.  over  the  P.F.I.  for  a  number  of  linear  proRramming 
problems  and  Bray  ton,  Gustavson  and  Willoughby  were  able  to  prove  this 
for  general  sparse  matrices. 

The  elimination  form  of  inverse  is  now  known  to  be  implemented  in  at  least 
two  commercial  linear  programming  codes;  Standard  Oil's  M5  code  and 
Scientific  Control  Systems'  UMPIRE  code  (see  Beale  C3}).  Other  codes  are 
in  the  course  of  being  modified  and  may  have  been  already.  Some  aspects 
of  the  P.F.I.  and  E.F.I.  inversion  techniques  will  be  examined  in  the  next 
section . 

The  success  of  the  Gaussian  or  Elimination  form  of  inverse  naturally  leads 
one  to  ask  whether  the  resulting  triangular  factors  can  be  used  to  advantage 
in  updating  the  inverse  in  the  iterations  following  a  re-inversion.  The 
answer  would  certainly  appear  to  be  yes,  Dantzlg  C8J  advocated  updating 
the  triangular  factors  of  the  basis  for  his  staircase  algorithm  since  this 
preserves  the  staircase  form  of  the  factors.  It  seems  extremely  likely 
that  in  this  special  case  updating  the  factors  rather  than  using  the  ordinary 
product  form  of  updating  would  lead  to  an  increase  in  efficiency.  Bennett 
and  Green  pursued  this  technique  for  general  sparse  matrices  and  Bartels  Cl) 
and  Bartels  and  Golub  C23  have  advocated  a  form  of  Dantzig's  method  on  the 
grounds  of  numerical  stability,  without  however  commenting  on  the  sparsity 
implications.  Brayton  et  al  C5)  have  proposed  two  other  updating  techniques 
although  their  concern  is  not  primarily  with  efficiency  in  terms  of  the 
simplex  method. 

The  three  new  proposals  will  be  reviewed  and  compared  in  Section  3. together 
with  their  implications  for  the  simplex  method.  One  of  these  methods  is  then 
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selected  for  some  computational  experiments,  the  results  of  which  are  presented 
in  Section  4.  Section  5  describes  how  this  method  might  be  Incorporated 
efficiently  into  a  production  code. 
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2. 


TllE  PRODUCT  AND  ELIMINATION  FORMS  OF  INVERSE 


P.F. I. 

A  popular  method  of  carrying  out  re-inversion  (see  Orchard-Hays  M)  is 
to  order  the  rows  (implicitly)  and  the  columns  (explicitly)  of  B  and 
partition  as  follows:- 


B  = 


B 

E 


A 


2 


where  A  and  A  are  lower  triangular.  The  submatrix  B  is  sometimes 

^  A 

referred  to  as  the  "bump"  and  those  vectors  with  elements  in  A  are  "above 
the  bump",  those  in  A g  "below  the  bump".  It  seems  to  be  characteristic  of 
most  linear  progtams  that  the  submatrices  A .  and  particularly  A.  are  quite 

X  A 

substantial  portions  of  B. 


The  above  partitioned  form  of  B  may  be  factorized  thus 


P  m 

I 

"l 

*1 

B  s 

A  I 

B 

I 

1 

D  I 

I 

E  I 

^2. 

Since  A,  and  A„  are  triangular  the  calculation  of  the  product  form  of 

X  « 

inverse  is  non-trivial  only  for  the  submatrix  B.  Various  "merit"  schemes 
have  been  proposed  for  selecting  the  order  of  pivots  in  B  but  they  will  not 
be  discussed  here  (see  e.g.  Tewarson  Cisj). 

Using  the  notation  of  reference  5  we  may  then  express  B  as  a  product  of 
elementary  transformations 
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D  =  Tj,  ... 


or  equivalently 


... 


where  T,  is  of  the  form  (we  assume  B  is  m  x  m) 
k 


Ik 


ft 

1  t 

^  I 


kk 

I  » 

j  'l 

4 


mk 


'  I  ^  (‘k  -  »k 


where  e  is  the  k^^  unit  m  vector,  e*  its  transpose,  and 

K  Iv 


0,-1  0,-1  V. 

K  =  T..  1  •••  T,  b, 

k  k-1  21k 


th 

where  b  is  the  k  column  of  B. 
k 


Note  however  that  two  transformations  are  made  for  each  column  in  the  bump  , 


hence  the  elements  of  the  for  such  columns  are  zero  on  the  rows  correi5ponding 


to  A,  and  in  addition  a  second  transformation  with  a  unit  pivot  and  tho  k 
column  of  E  on  the  rows  of  A2/IS  also  made.  Note  also  that 

=  I  -  tj^  (tj^  -  and  that  it  is  essentially  necessary  to  keep 

only  the  vector  t  to  store  T  ^and  t  and  t  J  to  store  T  (see  C53). 

K  K  K  KK  K 


th 
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K.F.  I. 


Ttie  abovt'  partitioninj;  schoino  may  be  used  with  only  a  slight  modi I'ioa ti on 
in  the  elimination  focm  of  inverse.  Uy  reai  ran^tinK  the  rows  and  e.oiuinns 
of  the  submatricos  ^2’  ^  repartition  B  to  give  (seo  Bealo  C^) 


D  = 


D 

A 


E 

n/ 

B 


/V 

where  A  2  is  now  upper  triangular.  This  procedure  obviously  facilitates 
factorization  into  triangular  factors  by  Gaussian  elimination  since  again 

A* 

only  the  factorization  of  B  is  non-trivial. 


Suppose  B  is  now  factorized  into  LU  where  L  is  lower  and  U  is  upper  triangular 
^  th 

Let  and  u  be  the  k  columns  of  L  and  U,  then  defining 

K  K 

Lr  =  I  +  ®k 

U  =  I  +  (u  -  e  )  e' 

k  k  k  k 


wo  may  factorize 


L  =  L  L_  • • •  L 
12  m 

U  =  U  U  ,  ...  U, 
ra  ra-1  1 


Hence 


0-1  .,-1  1-1  T-1 

B  =U  ...U  L  ...L 

1  mm  1 


There  is  some  latitude  in  the  choice  of  values  of  diagonal  elements  of  L  and  U 


The  most  efficient  choice  would  appear  to  be  setting  u^^^^  to  1  for  k  in 


pa 


rtitions  A.,  ^  and  to  1  in  partition  A..  If  this  is  done  the  number 

1  kk  2 

of  transformations  (disregarding  unit  transformations)  is  exactly  the  same  as 


for  the  P.F.I.  since  two  transformations  are  calculated  only  for  the  columns 


in  B. 
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The  order  of  choice  of  pivot  elements  in  B  may  be  decided  here,  as  for 


the  P.F.I.,  by  merit  schemes.  These  are  exhaustively  discussed  by 
Dantzig  et  al  c^. 


:j.  methods  of  updating  the  inverse 


In  this  st'Ctiori  we  briofly  roviow,  without  rigorous  proof,  and  compare  four 
methods  of  updating  the  inverse.  All  except  the  first  assume  the  elimination 
form  of  inverse  discussed  in  the  previous  section.  With  the  exception  of 
method  II  all  are  described  at  greater  length  by  Drayton  et  al  Csj. 

Method  1  C73 

This  is  the  standard  procedure  used  in  the  product  form  of  the  simplex  method. 
Suppose  column  b  of  the  basis  is  to  be  replaced  by  b  .  Let  B  be  the  original 

K  K 

basis  (factorized  in  either  E.F.I.  or  P.F.I.  form)  and  let  D  be  the  new  basis. 


Let 

then 

B-' 

where 

Hence  if 

B-' 

then 

'^/+i  •••  ^1 

Further  changes  are  made  by  repeating  the  above  process. 


Method  II  _L1 

This  is  Bartels's  version  of  Dantzig's  algorithm  C83.  Suppose  we  have 


B  =  LU 
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and  we  again  change  to  Then  L  ^  B  will  bo  identical  to  U  except 
in  column  k.  Permute  columns  k>l,...,m  one  place  to  the  left  and  place 
column  k  in  position  m  to  give  a  Hessenberg  matrix 


H  = 


=  l“^  B  Q 


(where  Q  is  a  permutation  matrix). 

The  elements  below  the  diagonal  in  columns  k  to  m-1  of  II  must  now  be 
eliminated  by  a  set  of  elementary  transformations.  An  important  feature  of 
Bartels*  method  is  that  in  carrying  out  the  elimination  rows  may  be  inter¬ 
changed  or  permuted  to  place  the  maximum  of  the  elements  h..  ,  h.  .  on 
the  diagonal  (Wilkinson's  "partial  pivoting"  strategy  Ciej)  to  ensure 
numerical  stability. 

The  result  is  a  new  upper  triangular  matrix 


U  =  G  ,  P  ,  ...  G,  P,  ,  G,  P,  H 
m-1  m-1  k+1  k+1  k  k 


where  the  P^  are  either  unit  matrices  or  the  permutation  matrices  exchanging 
rows  ^  and  /+1,  and  G^  are  the  elementary  transformation  matrices  eliminating 
the  element  below  the  diagonal  in  column  /  . 

Hence  b”^  =  Q  u"^  G  ,  P  ,  . . .  G,  P,  l"^ 

m-1  m-1  k  k 
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Dantzig  fs)  and  Bennett  and  Green  C4)  have  shown  that  the  "near  commutativity" 
of  the  matrices  G|  may  be  used  to  calculate  a  new  lower  triangular  matrix  L 
such  that 


l"^  =  G  ,  P  ,  G,  P,  l"^ 

m-1  m-1  k  k 


Thus  maintaining  a  true  triangular  decomposition.  This  is  not  however  an 
essential  feature  of  the  method. 

Of  course  the  permutation  matrices  ,  Q  need  not  appear  explicitly.  All 
that  Is  necessary  is  to  record  the  order  of  pivoting.  The  process  Is 
repeated  for  the  next  Iteration  using  the  new  U  and  (implicitly)  the  new  L. 


Method  III  J[5l 


Neither  this  method  nor  the  next  maintain  a  true  triangular  decomposition, 
however,  both  require  initially  the  elimination  form  of  Inverse. 

Suppose  we  have  B  =  LU  or  equivalently 


b"^  =  ...  u'^  l"^  ...  l7^ 

1  mm  1 


and  again  change  column  b  to  b  and  B  to  B. 

K  K 


The  method  proceeds  as  follows 


Let 


V  ::  L  ^  b 
k  k 


then  L  ^  B  is  identical  to  U  except  for  column  k  which  is  now  v  instead  of 

K 
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To  reduce  L  D  back  to  upper  Iriani'iilur  form  Ihi'  firsl  slop  ciirr iod  oul 
is  to  reduce  row  k  to  zero  In  columns  k-t-l,...,m. 


The  row  transformation  is  then  the  elementary  matrix 


W.  =  I  +  w ' 
k  k  k 


and  the  desired  elimination  Is  performed  by  forming  W  ^  L  ^  B,  since  the 

K 

columns  k-fl,...,m  of  L  ^  B  are  Identical  to  U,  W  ^  U  b  (I  -  e  w')  U  and 

iC  K  iC  ^ 

w'  U  is  simply  the  row  (0,...,0,u  ).  Furthermore  the  first  k 

K  K^K^X  KTH 

columns  of  W,  are  unit  vectors  and  hence  the  first  k  columns  of  L  ^  B  are 
k 

unchanged. 

Now  let  t  =  W  ^  L  ^  b  and  form  the  elementary  column  transformation 

K  K  K 


\  ^  ®k^  ®k 


then  clearly  forming  T^^  L  ^  B  reduces  column  k  to  the  unit  vector 

and  columns  k+l,...,m  are  unchanged  since  they  have  a  zero  on  the  pivot  row. 

(k) 

The  result  then  is  a  new  upper  triangular  matrix,  say  U  ,  obtained  from  U 


bll 

by  replacing  the  k  row  and  column  by  the  unit  vector  Oj^. 
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1.0. 


k  k 


or 


5-1  _  ^(^,-1  -1  -1  ^-1 
k  k 


(k) 

Using  product  form  notation^ if  we  lot  denote  the  elementary  matrix 

with  u  set  to  zero  we  may  re-write  B  ^  as 

K|f 


b"^  =  u"^  . . .  u"^ 

1  k-1  k+1 


y(k)-l  -1  -1  ^-1 

m  k  k  m 


Note  that  the  row  transformation  can  be  represented  (if  desired)  by  a 

product  of  two-element  elementary  column  transformations.  Also  note  that 

although  we  have  introduced  two  new  transformations  that  if  U  was  previously 

(k) 

in  product  form,  U  =  U^  ...  U^,  the  new  form  U  is  obtained  by  simply 

deleting  U  and  a  single  element  at  most  from  each  U,  ,  ...  U  for  a  net 
k  k+1  m 

gain  of  one  transformation  in  general  (unless  U  happens  to  be  the  unit  matrix). 

K 

To  carry  the  process  on  one  simply  treats  T  ^  W  ^  L  ^  as  the  lower  triangular 

Jv  )v 

factor  as  before,  although  in  fact  if  multiplied  out  the  resulting  matrix 
would  not  in  general  be  triangular. 


Method  IV  JSl 


This  method,  which  also  assumes  the  elimination  form  of  inverse,  proceeds  by 

inserting  matrices  T,  into  the  list  U  ...  U,  and  sometimes  deleting  a  U.  . 

k  ml  k 


To  change  b  to  b  we  form  the  vector 

K  K 


=  u 


-1 

'k+1 


m  m 


L-^  5 

1  k 


and  )(.‘place  U^^  by  the  transformation 


"^k  =  ^  ^  ®k 
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Then 


b“^  =  u7^  ...  u"^.  t’^  u“^,  ...  u"^  l”^  ...  l"^ 

1  k-lkk+1  mm  1 


Continuation  of  this  process  is  rather  more  cumpllcited  than  In  the  preceding 

methods.  Suppose  now  we  have  changed  b  to  b  and  wish  to  change  b.  to  b.  . 

k  K  £,  JL 

There  are  two  cases: - 


1.  If  £  ^  ^  then  let 


t  =  U, 


-1 


-1 


/+!  • 


.  U  "  t”^  U~^ 

k-1  k  k+1 


11“!  i-l  K 

U  Jj  •••Li  D 

mm  1 


form  T£  s  I  +  (t^  -  e^)  e^ 


and 


replace  by  T^  as  before. 


2.  If  £  >  k  form 


t|i  =  u’l,  ...  u"^  l“^  ...  l"^  b 
I  k  k+1  mm  1 


and  the  transformation  Tj^  as  above.  However,  is  not  replaced 


-1 


,-l 


and  Tj  is  placed  directly  before  T  in  the  inverse  transformation 


list. 


i.e.  b"^  = 


u"^  T~^  t”^  U~^ 

k+1  /  k  k+1 


u"^  ...  u'^  l"^  ...  l"^ 

1  m  1 


Continuing  further,  if  columns  kj^,...,kp  have  been  changed  and  column^  is 
to  be  altered,  let  k  =  min  (k^^  ...  k^)  and  carry  out  the  above  procedure. 


Having  completed  our  summary  of  methods  for  updating  a  basis  inverse  or 
substitute  inverse  we  must  consider  their  relative  merits  in  the  context  of 
the  simplex  method  in  terms  of  preserving  sparsity,  work  required  per 
iteration,  and  storage  and  data  handling  considerations. 
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J 


The  groat  advantage  of  Method  I  of  course  is  its  simplicity  and  the  minimum 
of  computational  stops  and  housekeeping  involved.  The  essential  step  in 
changing  the  basis  in  the  simplex  method,  given  the  constraints 


Ax  =  b 


and  the  now  column  a  to  enter  the  basis,  is  the  calculation  of  the  updated 


column 


OC  =  a 


and  the  ratio  test  with  the  updated  right-hand  side, 

^  b, 

to  find  the  pivot  row  k  from  (assuming  ^  0) 


e  s 


/^i 

min  '  — 
i 

^>0 


The  vector  ,  which  we  cannot  avoid  calculating,  is  precisely  the  vector 
t^^j^  used  in  Method  I  to  update  the  basis  inverse  (and  the  right-hand  side). 

Methods  II  to  IV  all  suffer  from  the  disadvantage  that  they  require  not  only 

the  vectored  for  determining  the  pivot  row  but  also  another  vector  t^^  which 

is  only  a  partially  updated  form  of  a..  Methods  II  and  III  can  more  or  less 

d 

overcome  this  disadvantage  by  putting  the  vector  L  ^  a  .  in  temporary  storage 

O 

in  the  course  of  calculating  oC  ,  at  the  expense  of  either  using  more  of  the 
computer's  core  storage  (thus  effectively  limiting  the  number  of  vectors 
which  can  be  saved  for  multiple  pricing)  or  of  using  backing  store.  This 
burden  is,  however,  comparatively  light  compared  to  the  immediate  and  acute 
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difficulty  encountered  in  Method  IV  where  the  vector  to  be  calculated 
for  updating  purposes  cannot  in  general  be  known  until  after  the  ratio  test 
operation  to  determine  the  pivot  row.  This  implies  that  t^^  must  therefore 
be  calculated  from  scratch  or  obtained  by  backward  transformations  on  the 
vector  cC  • 

Now  let  us  consider  the  implications  of  Methods  II  to  IV  in  terms  of  work 
per  iteration  expended  in  updating  the  basis  inverse  and  of  the  necessary 
data  manipulation.  Current  packing  methods  of  storing  sparse  data,  parti¬ 
cularly  strings  of  transformations  (see  C6},  tl5I  )  make  it  rather  difficult 
to  modify  an  existing  product  form  of  L  and  U.  Replacement  of  U  by  a  T 

K  K 

with  more  non-zeros  in  Method  IV,  for  exanple,  necessitates  either  a  whole¬ 
sale  shifting  of  elements  of  U,  maintenance  of  a  second  transformation  file 
to  be  inter-woven  with  the  first,  copying  out  a  complete  new  basis  file  or 
initial  and  wasteful  leaving  of  space  for  such  eventualities.  A  similar, 
in  fact  worse,  situation  occurs  in  Method  II  where  on  the  average  half  of 
the  transformations  must  be  modified  at  each  iteration,  increasing  the 
number  of  non-zeros  to  be  stored.  Explicit  storage  of  U  can  only  be  an 
answer  for  very  small  problems.  Some  kind  of  two-file  system,  each  taking 
the  new  inverse  in  turn,  would  appear  to  be  necessary  in  general,  involving 
a  great  deal  of  read/write  activity.  Method  III  on  the  other  hand  requires 
only  the  elimination  of  previously  non-zero  elements  in  U  -  a  comparatively 
easier  task.  This  point  will  be  returned  to  in  Section  5. 

The  amount  of  arithmetic  work  per  iteration  would  appear  to  be  greatest  in 
Methods  II  and  III.  The  bulk  of  the  work  in  Method  III  comes  from  applying 
the  transformations  to  H  to  reduce  it  to  triangular  form.  On  the 
average  we  expect  to  have  to  deal  with  ra/2  columns  on  the  right-hand  side  of 
the  matrix.  This  is  a  considerable  amount  of  computation  and  becomes  pro- 
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grossivoly  worse  as  U  fills  in  with  each  iteration.  For  Method  III  the 
main  effort  is  in  computing  the  row  vector  w'.  It  will  bo  shown  in 

K 

Section  5  that  this  calculation  may  be  performed  concurrently  with  the 
backward  transformation  phase  of  the  next  simplex  iteration.  The  most 
important  point,  however,  is  that  U  becomes  progressively  less  and  loss 
dense  in  this  method  and  the  calculation  progressively  easier.  There  is 
admittedly  the  extra  complication  in  Method  III  that  row  transformations 
are  called  for,  but  this  is  more  of  an  inconveneience  than  a  problem  and 
as  already  pointed  out  could  be  circumvented. 

Finally  we  consider  the  most  important  point  -  preservation  of  sparsity. 

In  their  discussion  Brayton  et  al  concluded  that  Method  III  was  superior 
to  Method  IV  in  this  respect  since  both  a  row  and  column  of  U  are  eliminated 
in  Exchange  for  a  vector  w',  which  only  has  elements  in  positions  k+l,...,m, 

K 

and  a  transformation  vector  t^^  which  is  the  incoming  column  multiplied  only 
by  L  On  the  other  hand  if  in  Method  IV  some  columns  with  very  small 

pivot  row  index  k  is  introduced  then  virtually  every  change  thereafter  will 
be  of  type  2  and  we  cannot  expect  much  improvement  over  the  ordinary  product 
form  method. 

In  comparing  Methods  II  and  III  it  again  seems  almost  certain  that  Method  III 
will  win  out.  The  number  of  new  elements  added  in  the  new  transformations 
would  appear  to  be  much  the  same,  but  as  already  mentioned  in  Method  II  U 
becomes  more  and  more  dense  while  in  Method  III  the  density  of  U  declines. 

In  this  connection  it  should  be  pointed  out  that  the  columns  of  U  most 
frequently  operated  upon  are  the  right-most.  These,  however,  are  just  those 
columns  which  will  be  initially  the  most  dense  after  inversion j since  inversion 
schemes  generally  postpone  the  densest  columns  until  last. 
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In  the  light  of  these  considerations  it  appears  that  Method  III  is  the 
most  promising  for  incorporation  into  the  simplex  method,  both  on  the 
grounds  of  work  per  updating  and  preservation  of  sparsity.  However, 
there  is  certainly  some  increase  in  updating  work  per  iteration  over  the 
standard  product  form  (Method  I)  and  a  significant  improvement  in  growth 
of  non-zeros  in  the  basis  representation  is  required  to  justify  its  use. 
Some  experimental  results  are  presented  in  the  next  section. 

Although  we  have  chosen  Method  III  ao  the  best  for  our  purposes  this  is 
not  to  say  that  the  others  may  not  be  si^erior  in  other  circumstances. 
Bartels's  method  (II)  is  of  considerable  importance  in  achieving  highly 
accurate  solutions  to  small  dense  problems  (see  Bartels  and  Golub  C23) 
while  Dantzlg's  original  version  of  Method  II  for  specially  structured 
matrices  is  as  yet  untried. 
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1. 


C0MPUTAT10NA.L  RESULTS 


Some  compututional  oxperimonts  have  been  carried  out  to  determine  whothor 
Method  III  gives  a  sufficient  improvement  over  Method  I  in  terms  of  growth 
of  non-zeros  to  warrant  further  investigation. 

To  carry  out  these  experiments  two  linear  programming  codes  were  written 
and  run  on  the  Stanford  University  IBM  360/67.  Both  are  all-in-core 
FORTRAN  codes  using  the  elimination  form  of  inverse.  Neither  employs 
multiple  pricing.  The  first  code  updates  the  Inverse  in  standard  product 
form  fashion  (Method  I),  the  second  is  a  modified  version  using  Method  III 
to  update  the  Inverse.  No  attempt  has  been  made  to  evaluate  timings, 
dependent  as  they  arc  on  machine  used  and  programming  technique.  Only  the 
number  of  non-zero  elements  in  the  inverse  representation  is  considered. 

Five  problems  ranging  from  small  to  moderate  have  been  used.  These  problems 
arc  two  of  the  well-known  SHARE  standard  test  problems,  26  and  IB,  (see  Smith 
and  Orchard-Hays  1143  )  and  3  other  test  problems  supplied  by  Joel  Cord  of  IBM. 
The  relevant  statistics  on  these  problems  are  given  in  Table  1. 

Only  the  smallest  problem  was  run  to  optimality,  the  remainder  being  cut  off 
after  200,  250,  300  and  350  Iterations.  To  compare  the  two  methods  the  set 
of  problems  was  run  on  both  codes  with  a  fixed  inversion  frequency  -  40 
iterations  for  the  SHAKE  problems  and  50  for  the  others.  The  initial  com¬ 
parison  is  arrived  at  by  counting,  at  Intervals  of  10  iterations,  the  number 
of  new  non-zeros  added  to  the  inverse  since  the  last  re-inversion.  To 
further  remove  constant  factors  only  non-pivot  elements  are  counted.  (Both 
methods  give  a  net  addition  of  one  pivot  element  per  iteration  in  general.) 
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All  problems  were  started  from  a  slack  basis  and  hence  the  iterations  up  to 
the  first  re-inversion  are  Ignored.  The  number  of  new  non-pivot  non-zero 
elements  after  multiples  of  10  iterations  from  re-jnverslon  is  then  averaged 
for  the  remainder  of  the  run.  These  figures  are  given  in  Table  2. 

The  raw  data  appears  in  more  digestible  form  in  Table  3,  where  Method  I  has 
been  used  as  a  base  and  the  fraction  of  new  elements  produced  by  Method  III 
Is  given.  The  improvement  is  obviously  very  satisfactory,  with  the  possible 
exception  of  problem  2.  As  might  be  expected  the  Improvement  tails  off  as 
the  number  of  iterations  Increases,  particubrly  for  the  smaller,  denser 
problems. 

To  obtain  a  better  estimate  of  the  effect  of  Method  III  on  the  simplex  algorithm 
we  should  also  consider  the  total  number  of  non-pivot  non-zeros  in  the  inverse. 
To  this  end  we  list  the  Inversion  statistics  in  Table  4  and  compute  the 
fraction  of  total  non-pivot  non-zeros  produced  by  Method  III,  again  using 
Method  I  as  a  base.  The  inversion  statistics  do  not  Include  the  original 
all-slack  basis  and  are  of  some  interest  in  themselves.  The  pivot  choosing 
procedure  in  decomposing  B  (see  Section  2)  is  essentially  algorithm  7  in 
Dantzlg  et  al  Csl  . 

The  modified  data  for  Method  III  appear  in  Table  5, 

Although  the  problems  are  comparatively  small  and  the  sample  is  also  small 
the  results  in  Table  5  are  quite  encouraging.  It  appears  that  for  reasonably 
sparse  problems  with  200-300  constraints  a  reduction  of  20-30%  in  the  number  of 
non-zeros  in  the  substitute  inverse  may  be  expected  using  Method  III,  thus 
considerably  reducing  the  time  spent  in  forward  and  backward  transformation 
in  the  simplex  method.  The  results  for  problem  4  show  that  the  reduction 
may  be  dramatic. 
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IMI'miKNTATlON  OF  METHOD  111 


Tho  rL'sults  i)f  tho  pri'C'Hlin^;  si'Ction  suggest  that  further  experiments  using 
a  modified  production  or  commercial  code  would  be  worthwhile.  It  would  also 
appear  that  the  modifications  could  be  made  and  the  extra  updating  work  per 
iteration  accomplished  at  comparatively  little  cost. 

Tho  first  essential  is  clearly  to  have  a  code  which  employs  the  elimination 
form  of  inverse.  Given  such  a  code  it  should  bo  a  simple  matter  to  adjust 
the  storage  and  buffering  activities  such  that  the  T  ^  W  ^  transforms  can  be 

IV  IV 

adjoined  to  the  left  of  the  product  form  L  ^  .  L,  ^  of  L 

m  1 


th 

The  reduction  of  tho  k  row  and  column  of  U  to  the  unit  vector  e  can 

k 

probably  be  accomplished  satisfactorily  by  simply  setting  indicator  bits. 

Usually  we  will  have  U  ^  represented  as  U,^  U  ^  ...  U  or  rather  the  packed 

14  m 

vectors  u^  corresponding  to  these  elementary  matrices.  The  reduction  may 
then  bo  accomplished  by  tagging  with  a  bit  to  indicate  it  is  to  be  skipped 
over.  Furthermore  using,  say,  the  blocked  index  scheme  (see  Smith  D3)  for 
storing  the  vectors  u^  tho  row  index  k,  if  it  appears,  may  be  tagged  to 
indicate  the  element  is  zero  for  k+l,...,m.  Failing  this  the  u^^^  elements 
could  bo  physically  replaced  by  zeros  if  necessary. 


This  scheme  is  of  course  inefficient  in  the  sense  that  in  an  out-of-core 
system  progressively  more  and  more  valueless  data  must  be  buffered  in  and  out 
of  core.  On  the  other  hand  if  the  system  is  already  compute  bound'. rather 
than  I/O  bound  this  is  no  restriction,  and  even  if  this  is  not  the  case  the 
very  frequent  re-inversions  carried  out  for  large  problems  should  prevent 
this  inefficiency  from  reaching  really  noticeable  proportions. 


Finally  let  us  consider  the  process  of  computing  w'.  For  simplicity  let 

Iv 

us  assume  we  are  carrying  out  the  first  Iteration  after  a  re-lnvcrslon. 

We  then  have 


Now  the  calculation  of 


-1 


and  also  the  modification  of  U  to  may  be  carried  out  in  a  single  pass 

from  left  to  right  through  the  transformations  U,^  ...  U  To  see  this 

1  ID 

notice  that  in  computing  say 


. V  '‘I 


-1 


for  £  >  k 


only  the  elements  play  part  in  the  arithmetic  since  U  and 

are  upper  triangular.  The  elements  z.  z  are  unchanged.  We 

c  L+1  m 

f  k) 

therefore  arrange  the  computation  of  w^^  and  the  modification  of  U  to  U 
as  follows:- 


Initial  Step 

Create  a  zero  row  vector  and  skip  transformations  U  ^  u”^  ,  Tag  u”^. 

1  k-1  k 

Let  £  =  k+1. 

General  Step 

For  transformation  ,  if  is  non-zero  extract  it  and  add  it  in 

position  £  of  the  row  vector,  marking  the  element  as  now  being  zero  in 
Post-multiply  the  row  vector  by  Proceed  for  £-  k+l,...,m. 
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\  vory  uilvuntagi'ous  result  of  this  procedure  is  that  the  updating  of  the 
inverse  can  be  carried  out  c(>ncurrently  with  the  backward  transformation 
of  the  following  simplex  iteration.  In  calculating  the  pricing  vector  for 
the  next  iterai ion  we  must  compute 

tr'=  c'  l‘^ 

k  k 

(k) 

for  some  row  vector  c'.  Now  the  first  k-1  columns  of  U  are  identical 
with  those  of  U.  The  modification  of  the  remaining  columns  of  U  may  be 
carried  out  in  the  same  pass  through  the  transformation  file  as  that  for 
the  calculation  of  Tf '  and  furthermore  only  one  unpacking  of  the  non-zeros 

from  each  u.  is  necessary.  Having  obtained  w'  the  new  transformations 

C  k 

may  be  attached  to  L  ^  and  the  remainder  of  the  calculation  of  if' 
can  proceed.  Looking  to  the  future,  this  procedure  is  admirably  suited 
to  parallel  processing  facilities. 
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6. 


CONCLUSION 


Linear  progranmilng  Inversion  routines  have  now  reached  a  very  high  level  of 
sophistication  and  efficiency,  both  in  terms  of  speed  and  sparseness  of  the 
resulting  inverse.  Although  there  is  doubtless  still  some  room  for  improve¬ 
ment  it  seems  likely  that  further  improvements  in  efficiency  must  be  sought 
in  other  areas.  The  attempt  to  maintain  a  sparse  inverse  is  one  such  area 
and  the  results  of  Section  4  give  considerable  grounds  for  optimism.  What 
is  needed  now  is  full  scale  experimentation  with  some  of  the  methods  reviewed 
here  on  large  problems  of  1000  and  more  constraints  to  confirm  and  extend 


these  results. 
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